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ABSTRACT 



The application of digital computers for the simulation of physical 
systems has become widespread. This paper describes a program designed 
to measure the frequency response of the simulation model of a linear 



system. 
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I. INTRODUCTION 



The application of digital computers for the simulation of physical 
systems has become widespread. Computer programs have been developed 
that simulate mechanical, electrical, and electro-mechanical systems, 
circuits, chemical processes, biomedical problems, traffic control, and 
so on . 

One of the most important groups of programs for the simulation of 
continuous systems is that which employs digital-analog simulators; 
i.e., programs which simulate the elements and organization of the 
analog computer. 

The purpose of the investigation presented in this paper is to 
develop a technique -- a digital computer program -- for measuring the 
frequency response of the simulation model of a system. While such a 
program should be applicable to any digital simulation language, the 
program description presented in this paper utilizes International 
Business Machines Company's DSL/360 Digital Simulation Language. 

DSL/360, a System/360 FORTRAN IV program for the digital simulation 
of continuous system dynamics, employs the building-block approach of 
digital-analog simulators while providing the power of logical and 
algebraic equation notation. 

DSL/360 provides a basic set of function blocks from which a 
physical system may be modeled: integrators, limiters, pulse generators, 

function generators, and so on. In addition, FORTRAN library functions 
and functions from the Scientific Subroutine Package (SSP) may be 
utilized. In the event none of these satisfies the user's needs, the 
user may provide his own function blocks or subroutines. 
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A digital computer program to determine the steady-state response of 
a system to a sinusoidal input is a useful adjunct to a digital simulation 



language. With 
be detected and 
Reference 1 



it, non-linearities in an ostensibly linear system may 
eva luated . 

describes in detail the use and operation of DSL/360. 
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II. DISCUSSION 



This chapter contains a discussion of the response of a linear 
system to a sinusoidal excitation and the basic concepts utilized in 
the development of a digital computer program to determine such response. 

A. SYSTEM RESPONSE TO A SINUSOIDAL INPUT 

When a system is linear, its response to a sinusoidal input is a 
sine wave of the same frequency (the higher harmonics are negligible). 

The magnitude ratio and phase of the response depend on the forcing 
frequency but not on the input magnitude [Ref. 2]. The condition that 
magnitude ratio and phase be independent of the input amplitude is a 
condition for the linearity of the system. 

The steady-state frequency response of a stable, linear system to 
a sinusoidal input can be determined analytically from the system 
transfer function [Ref. 3]. The response to an input A sin cut is 
given by 

y = a| P(jo)) | sin(o)t + §) 

where | P(j<x>) | is the magnitude of P(jo)), $ is the argument of P(jo)), 
and the complex number P(jo)) is determined from the system transfer 
function P(s) by replacing the s with jo). The system output has the 
same frequency as the input and can be obtained by multiplying the 
input by P(jo)) and shifting the phase angle of the input by the argu- 
ment of P(jo)). | P(jao) | and $ for all (jd constitute the system frequency 

response, where P(jo)) is the gain of the system for sinusoidal inputs 
with frequency a). 
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B. RESONANCE IN A SYSTEM AND STABILITY 



The maximum value of the magnitude of the closed-loop frequency 
response of a system is a measure of the stability of the system 
[Ref. 3], The frequency at which this maximum occurs is the resonant 
frequency of the system. 

It often occurs in reality that a linear system’s response will 
contain non-linearities under certain conditions. As an example, a 
second-order system may experience system gain non-linearity as a 
result of resonance if the excitation frequency is at or near the 
system’s natural frequency. Such non-linearity and any resultant 
effect on system stability are of interest to a design engineer, for 
instance, since the performance of the system with which he is con- 
cerned may be vitally affected. Being able to utilize a computer to 
determine system non-linearity during the design stage may result in 
savings of both time and money. 

C. DEVELOPMENT OF THE PROGRAM 

A program to determine the frequency response of a system to a 
sinusoidal input must first determine attainment of the steady-state 
condition upon application of the excitation and then calculate values 
of the steady-state response magnitude and phase. 

D. ASSUMPTIONS 

Several : assumptions are made regarding the system, its excitation, 
and its response. 

Assumption 1. The system is linear and stable. 

Assumption 2. The system may be described in DSL/360 format. 

Assumption 3. There are no initial conditions; i.e., all initial 
conditions are zero at time zero. 
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Assumption 4. The excitation is described by 

A sin(o)t 4* 9) 

where A is the peak amplitude, co is the angular frequency in radians 
per second, and 9 is the phase (input phase is zero throughout this 
paper) . 

Assumption 5. The steady-state output of the system is periodic 
with the same frequency as the input and may be described by 

B sin (cut + $) 

where B is the peak amplitude, a) is the angular frequency in radians 
per second, and $ is the phase of the output with respect to the input. 

E. DETERMINATION OF THE STEADY-STATE CONDITION 

In any digital simulation language, the independent variable is 
time. The computational process is an iterative one, performed at 
intervals of time specified by the user. Attainment of the steady- 
state condition is determined by a mechanical, va lue -c omparing process 
rather than by the solution of an equation, such as that discussed 
earlier . 

The steady-state condition is determined in the program discussed 
in this paper by permitting the system to respond to the sinusoidal 
input until certain steady-state criteria are met. Ten cycles of input 
are permitted to pass before attempting any steady-state testing, by- 
passing computations involving the more wildly fluctuating transients 
and, thereby, decreasing computation time. The figure "ten" is an 
arbitrary choice, based on the assumption that the output will have 
settled to a condition approaching the steady state after ten cycles 
of input. 
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The steady-state determination is accomplished by comparing the 
value of the output magnitude at each iteration point on one cycle 
with the value at the same point on a later cycle. When these values 
are reasonably close, steady state is considered to have been attained. 
To reduce the possibility of chance satisfaction of the steady-state 
criteria while the output is still in the transient or settling phase, 
comparison of output magnitude values is made in the program with not 
one but two later cycles, the third and the tenth. Again, these figures 
are arbitrary choices. 



F. DETERMINATION OF OUTPUT PEAK MAGNITUDE AND PHASE 

Once the steady-state testing commences, the program simultaneously 
computes the peak magnitude of the output and its phase with respect to 
the input. 

The peak magnitude of the output is nothing more than the highest 
steady-state value computed for the output magnitude. 

The phase is computed from the times at which an input cycle and 
its resultant output cycle cross their respective zero-value points. 

The equations used in the program are: 






2rr(t.-t ) 
1 o 






360(t.-t ) 

1 o 



where is the phase in radians, is the phase in degrees, t^ is the 
time at which the input cycle crosses its zero-value point, t Q is the 
time at which the resultant output cycle crosses its zero-value point, 
and T is the period of both the input and the output. 
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III. THE PROGRAM 



This chapter presents a description of the program evolved to 
determine the response of a linear system to a sinusoidal excitation. 

A sample program of a single-run job is shown in Figure 3.1. Reference 
to this program is made in the explanations in the paragraphs that 
follow. The explanations presuppose some familiarity with DSL/360 on 
the part of the reader. 

A. THE SYSTEM 

The system used in the sample program of Figure 3.1 is second-order, 
with a natural frequency of ten radians per second, but any linear 
system would have served as well. The system is described, from its 
closed-loop form, as having a forward-loop transfer function of 
100/s (s + 1) and a negative unity feedback- path transfer function. 

The system might also have been modeled in its open-loop form. 

The system modeled in this manner is shown in Appendix A. 

B. CONSTANTS, PARAMETERS, AND PROGRAM CONTROLS 

Constants, parameters, and program execution controls are shown in 
the sample program following the title and general system description 
(in DSL/360, asterisks in column 1 indicate a comment card). 

The one constant, PI, is self-explanatory. 

The parameters must be determined by the user and are: 

OMEGA - the input frequency, in radians per second 

A the input peak magnitude 

FGAIN - the feedback- path gain 
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TITLE LINFAR SYSTEM RESPONSE TO A SINE WAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: 100/SIS+l) 

* FEEDBACK-PATH TRANSFER FUNCTION: -1 
CCNST PI =3. 1415927 

PARAM OMFGA=2C. , A=l. ,FGAIN=1. 

CONTRL FINTIM=!CCOO. , DEL T=0« 015707964 
INTEG RKSFX 
INITIAL REGION 

PERI0D=2.*PI/0MEGA 
DELINT=3. ♦PERIOD 
DILINT=10.*PERI0D 
ERG=PI*(.5-DELT/ PERIOD) 

EPS1=1. E-03*A 
EPS2=0. 

EPS3=C. 

XIN=0. 
xout = o • 

XTDI FF=0. 

RPHASE=0. 

DPHAS E = 0. 

MAXOUT = 0. 

DEROUT = 1C. 

DERIVATIVE REGION 

I N= A*S I NF ( C. , OMEGA, 0, ) 

ERPOP=IN-FGA IN*OUT 

0UT=TRNFR(0. ,2. , IC , NUM, DEN, ERROP ) 

STORAG I C ( 2 ) , NUM ( 1 ) , D E N ( 3 ) 

TABLE ICC 1—2 ) =0. ,0.,NUM( 1) = 100. ,DEN<1-3)=1. ,1. ,0. 

DYNAMIC RFGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DELOUT = DEl AY ( 60 5, DEL I NT, OUT ) 

DI LOUT = DEL AY (2005, DILINT,OUT) 

IFITIME.LE.DELINT) GO TO 2 
I F ( ABS ( DEL OUT- OUT ).LE.EPS1) GO TO 1 
MAXOUT = 0. 

GC TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDF 

1 IFIOUT.GT. MAXOUT) MAXOUT=OUT 
DEROUT=OMEGA*MAXOUT*COS< ARG) 

EPS2=C.OC1*MAXOUT 
EPS3=0MEGA*MAX0UT*C0S( ERG) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 AXIN=CRCSS (TIME, IN,0. ) 

I F( A X IN • NE • 0* ) XI N = AX I N 
A XOUT=CRC SSITIME* OUT , 0 • ) 

IFIAXOUT.NE.O. ) XOUT=AXOUT 
IFIXIN.GT. XOUT ) GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE=2*P I*XTDIFF/PER I OD 
DPHASE=36C.*XT0IFF/PERI00 

3 ARG= OMEGA* TIME+R PHASE 
IFITIME.LE. DILINT) GO TO 4 

I F I ABS I DEL OUT-OUT ) . LE. EPS 2. AND. ABS I DI LOUT-OUT. , . 

) .LE.EPS2. AND. ABS I DEROUT. . . 

). LE.EPS3. AND. OUT. GT.O. .AND. MAXOUT. GT.O. ) CALL ENDRUN 

4 CONTINUE 
TERMINAL REGION 

PRINT OMEGA, MAXOUT ,RPHASE ,DPHASE 
CALL PRINT 

END 

STOP 



FIGURE 3.1 

SAMPLE PROGRAM FOR DETERMINING SYEADY-STATE RESPONSE TO A 

SINUSOIDAL EXCITATION 
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The program execution controls must also be determined by the user 



and are : 

FINTIM the maximum simulation value for the independent 

var iab le , t ime 

DELT the simulation interval; the unit of time for the 



integration routine to accomplish an integration step 

INTEG RKSFX - the fixed integration routine with an integration 
interval equal to one-half of the simulation interval, DELT 

FINTIM must be expressed numerically in DSL/360 and, in this program, 
is set arbitrarily at 10,000 seconds to ensure sufficient time for the 
response to attain the steady-state condition. Once steady-state is 
attained and the output peak magnitude and phase are calculated, the 
run is automatically terminated by a CALL ENDRUN statement. 

In this program, the execution control DELT is directly related to 
OMEGA, the input frequency. Since DSL/360 requires that DELT also be 
expressed as a numerical value, the DSL/360 user must calculate the 
value of DELT from the equation 

DELT = ^ x PERIOD = ; ; 2 seconds. 

N N x OMEGA 

The first iteration calculation is made at time zero, and an iteration 
calculation is made at discrete intervals thereafter, the length of the 
intervals depending on the integration scheme being used [Ref. 1]. 

Since the output phase is not known in advance, assurance that an 
iteration calculation will be made at or reasonably near the peak- 
value point of each output cycle is made by setting a high value for 
N. Since the basis for steady-state testing is the comparison of values 
at each iteration point on one output cycle with the values at the 
corresponding point on two other output cycles, N should be an integer. 

It has been determined empirically that the value of N should be an 
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integer no less than about twenty; this figure provides a trade-off 
between acceptable accuracy and reasonable computation time. A program 
for generating a table of DELT versus OMEGA, for N equal to twenty, is 
shown in Appendix B. 

C. THE INITIAL REGION 

The INITIAL REGION encompasses calculations, input and output 
operations, and initializations that must be made once only at the 
beginning of a run. The values and equations shown in the sample 
program apply to all systems and need not be changed by the user. 

All variables shown in the INITIAL REGION in the sample program 
are either self-explanatory or explained more appropriately elsewhere 
in this chapter. 

D. THE DERIVATIVE REGION 

The DERIVATIVE REGION encompasses those calculations involving 
integration and derivatives of the state variables being integrated. 

The basic interval in the independent variable, time, for each pass 
through this region is the calculation interval. In this program, the 
interval is determined by INTEG RKSFX and is one-half the simulation 
interval, DELT. 

In the program, the DERIVATIVE REGION contains equations describing 
the input and system transfer functions and input-output relationships. 
While the DSL/360 function block TRNFR is used in the sample program, 
representing the system in either its closed-loop or open-loop form, 
the system could also have been represented by a set of ordinary dif- 
ferential equations [Ref. 1]. The use of ordinary differential 
equations does not alter the remainder of the program in any way. 
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STORAG, a DSL/360 translator command, and TABLE, a DSL/360 data 
input, are both associated with the use of the function block TRNFR. 
SINEQ is the DSL/360 function block that models the input. 

E. THE DYNAMIC REGION 

The DYNAMIC REGION encompasses all time -dependent algebraic calcu- 
lations, other than derivative calculations for the integration routines 
in the DERIVATIVE REGION, plus all other operations which must be per- 
formed at each discrete value of time; i.e., all those operations which 
must be performed each iteration. It is in this region, in the sample 
program, that steady-state testing and response magnitude and phase 
calculations are accomplished. 

The variables in the DYNAMIC REGION of the sample program are: 

OUT the value of the output magnitude at the point on the output 

cycle at which the current iteration calculations are being made 

DELOUT - the value of OUT, to be saved for comparison with the value 
of OUT at the same point on the third cycle following 

DILOUT - the value of OUT, to be saved for comparison with the value 
of OUT at the same point on the tenth cycle following 

EPS1 the allowable difference between DELOUT and OUT (the differ- 

ence between the current value of the output magnitude and the value at 
the same point on the third cycle preceding) within which DELOUT and OUT 
are considered to be equal 

MAXOUT - the maximum value of OUT; the peak value of the output 
magnitude 

DEROUT - the value of the slope of the output at the iteration point 

EPS2 the allowable difference between DELOUT and OUT and between 

DILOUT and OUT within which DELOUT, DILOUT, and OUT are considered to 
be equal 

EPS3 the maximum allowable value within which the slope of the 

output, DEROUT, may be considered to be at the peak of an output cycle 

AXIN the time at which the input cycle crosses its zero-value 

point 
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XIN the value of AXIN saved until the input cycle next crosses 

its zero-value point 

AXOUT -- the time at which the output cycle crosses it zero-value 
point 

XOUT the value of AXOUT saved until the output cycle next crosses 

its zero-value point 

XTDIFF - the time elapsed between XIN and XOUT 

RPHASE - the value of the output phase in radians 

DPHASE - the value of the output phase in degrees 

ARG the value of the argument for the calculation of the output 

slope, DEROUT 

ERG the value of the argument for the calculation of EPS3 

EPS1 is relatively coarse and is the criterion to be met by the 
output magnitude before the program will permit steady-state testing 
during any given iteration. EPS2 is a refined criterion for determining 
whether or not the train of cycles is sufficiently similar for the out- 
put to be considered in the steady-state condition. The run is continued, 
however, even after the EPS2 criterion is satisfied, until the peak value 
of the output magnitude has been attained, determined, in part, by EPS3. 

The numbers 605 and 2005 appearing in the arguments of the equations 
for DELOUT and DILOUT represent the maximum number of sampled values of 
OUT stored in the delay intervals DELINT and DILINT, respectively. 
Reference 1 states that these numbers , which must be coded explicitly 
as numerical, integer constants, should be less than or equal to the 
delay interval divided by DELT. In practice, however, it was found 
that these numbers must be greater than the delay interval divided by 
DELT. For an N of twenty, then, in the case of DELOUT, this number is 

DELINT . 3 x PERIOD . 

DELT " 1/20 x PERIOD “ ° U ' 
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For DILOUT, the number is 



DILINT 10 x PERIOD s 
DELT 1/20 x PERIOD 

It can be seen, therefore, that the sample program will handle a pro- 
gram for any DELT where N is 204 or less. If more than this number of 
iterations per cycle is desired, the integer numbers to replace the 
numbers 605 and 2005 may be determined from the above equations. 

ERG is the argument, in radians, for the slope of the output at 
the point on the cycle most remote from the maximum point which can 
possibly be obtained for the DELT specified. Ideally, an iteration 
calculation will occur at the maximum, and the slope will be zero. 
However, because the process is an iterative one and the phase is 
shifted by an indeterminable amount, an iteration calculation may occur 
as far away, in time, as one-half of DELT; i.e., the maximum point may 
be exactly half-way between two iteration calculation points. To 
ensure that the slope criterion, EPS3, will be no greater than the slope 
at this point, ERG is determined as follows: 

Nr. of radians/iteration - 



2 x pi 



Nr. of iterations/cycle 
2 x pi 



PERIOD/DE LT 

_ 2 x pi x DELT 
PERIOD 



ERG 



» £i _ i r. 

2 2 L 

[o. 



2 x pi x DELT 
PERIOD 



= pi 



5 - 



DELT 

PERIODJ 



Phase is calculated in the program so as to produce the result 
always as a lag indication by XTDIFF, else both the true phase angle 
and its supplement would be calculated. For a phase angle from zero 
degrees to 180 degrees, the program is satisfactory. However, because 
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all attempts to limit phase computations only to times at which the 
input and output cycles crossed their respective zero-value points in 
an upward (increasing value) direction failed, phase lag angles from 
180 degrees to 360 degrees (phase lead) are calculated as supplements 
of the true phase angle. To determine whether the program has calcu- 
lated true phase or its supplement, the user may either sketch the 
asymptotes of the phase plot on a Bode diagram or, if the system 
description does not readily lend itself to this method, have the pro- 
gram provide a time graph of a few cycles each of the input and the out- 
put, from which phase lag or lead may be determined visually. 

F. THE TERMINAL REGION 

Entry into the TERMINAL REGION is made at the termination of a run 
for purposes of performing input and output operations, testing 
terminating conditions, processing results, changing parameter values 
and requesting a rerun, or terminating the job. In the sample program, 
the values of the variables and parameters desired in the printout are 
caused to be written in the TERMINAL REGION. 

If it is desired to make several runs in one job, each run at a 
different value of OMEGA, say, it is within this region that new 
values for the parameter OMEGA and the simulation interval DELT may 
be specified. Appendix A contains an example of such a multi-run job. 

A method for programming a multi-run job with different values for 
OMEGA for each run, calculated by the program, and for which DELT need 
not be changed for each run is discussed in Chapter IV. 
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IV. CONCLUSIONS 



This chapter presents a discussion of the results obtained, several 
limitations of the program as evolved, and recommendations for exten- 
sions of investigation. 

A. RESULTS 

Overall, the results of testing several linear systems were satis- 
factory. The response data provided by the program were consistently 
accurate, signifying that the steady-state testing scheme worked 
satisfactorily and consistently avoided chance satisfaction of the 
steady-state criteria during the settling phase. 

Programs and results for several of the systems tested are shown 
in Appendix A. 

B. LIMITATIONS 

There are several limitations on the use of the program evolved, 
some inherent in DSL/360 and some a function of the value of parameters 
se lected . 

Limitation 1. DELT must be calculated by the user and expressed in 
the program numerically. This could be eliminated if DELT could be 
expressed as a variable function of OMEGA, but this requires an altera- 
tion of DSL/360, beyond the scope of this paper. 

Basically, a new value for DELT must be calculated for each value 
of OMEGA used. However, the number of calculations for values of DELT 
may be reduced in a multi-run job by starting with a high value for 
OMEGA and a corresponding value for DELT (equal to, say, one -twentieth 
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of the period) and including in the program a scheme to halve the pre- 
ceding value of OMEGA each run. DELT need not be changed, then, since 
each time OMEGA is halved, the number of iterations per cycle is 
doubled, and the criteria for DELT are satisfied (N remains an integer 
no smaller than twenty, doubling each successive run; the program as 
shown in the sample will accept any value of N to a maximum of about 
204). Appendix A contains an example of this method. 

Limitation 2. Since the independent variable is time, only time 
plots may be obtained without modification of DSL/360. Such plots may 
be desired if the user wants to determine whether phase lags or leads. 

Limitation 3. Use of too small a value for OMEGA, as compared to 
the value for the natural frequency of the system, may result in no 
program output. This is caused by the fixed integration called for by 
INTEG RKSFX. If this card is removed from the deck, the integration 
interval automatically defaults to that called for by the variable- 
step integration routine RKS . However, DELT then varies with the 
integration interval [Ref. 1], and the conditions required for steady- 
state testing may not be met because integration points on different 
cycles may not correspond. It was attempted to begin a run with the 
variable-step routine RKS and, after some relatively short interval, 
to impose the fixed-step routine RKSFX. These attempts were unsuccess- 
ful, however, because no matter where the INTEG RKSFX card was located 
-- even within a N0S0RT region [Ref. 1] -- it was employed from the 
outset . 

C. RECOMMENDATIONS FOR EXTENSIONS OF INVESTIGATION 

There are several recommendations for logical extensions of this 
investigation. 
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Recommendation 1. Determine a method for calculating phase that 
will not result in computing the supplement of the phase angle, rather 
than the true phase angle, whenever the phase is leading (or, as 
written in the program, whenever the phase lags from 180 degrees to 
360 degrees) . 

Recommendation 2. Alter DSL/360 to permit the designation of DELT 
as a variable function of OMEGA rather than as a numerical constant. 

Recommendation 3. Alter DSL/360 to permit graphing a variable 
versus some other parameter than time; specifically, it would be 
desirable to graph both MAXOUT and the phase versus OMEGA. 

Recommendation 4. Apply the basic concept of this paper to the 
measurement of frequency response for a non-linear system. 
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APPENDIX A 



This appendix presents six examples of systems used to test the 
program evolved for determining the frequency response. For each 
example, the system description, the program, and the frequency response 
is given. 
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EXAMPLE 1 



The system used in this example is identical to the system used in 
the sample program of Figure 3.1. This example shows one method of 
obtaining the response to more than one frequency in one job. DELT 
must be calculated and specified within the program for each OMEGA. 

The program is shown in Figure A.l, and the frequency response is 
shown in Figure A. 2. 

For this example, the simulation time varied from 3.65 seconds to 
8.11 seconds for runs for nine different frequencies, or an average of 
6.05 seconds per run. The total execution time was 58.49 seconds. 
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TITLE LINFAF SYSTEM RESPONSE TO A SINE WAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: IOO/SIS+I) 

* FEEDBACK-PATH TRANSFER FUNCICN: -1 
CONST PI=3. 1415927 

PARAM OMEGA =6. »A=1 . ,FGAIN=1. 

CJNTRL F I NT I M= 1 OOOC. , DEL T = 0 . 0523 5936 
INTEG RKSFX 
INITIAL REGION 

PERIOD =2. *PI /OMFGA 
DEL IN T = 3. ♦PERIOD 
DIL INT=10.*PER IOD 
ERG = P I * ( . 5-DEL T /PER IOD) 

EPS1=1 .E-03*A 
EPS2=0. 

E°S3 = 0 . 

XIM=0. 

XOUT= G . 

XTDIF C =C. 

RPHASF=C . 

0PHASE=0. 

MAX0UT=0. 

DEROUT = 10. 

DERIVATIVE REGION 

IN = A* S INE ( 0. .OMEGA ,0. ) 

ERROR= I N-FGAIN*OUT 

Ol)T=TRN c P ( 0 • ,2. ,IC »NU M »DEN .tRRGR) 

STORAG IC( 2 ) ,NUM( 1 ) ,DEN( 3) 

TABLE IC( 1-2 )=0. ,0. ,NUM( 1 ) =100. , DE M ( 1 -3) = 1 . , 1 . ,0 . 

DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DELOU T = DELAY( 605 , DEL I NT.OUT ) 

DIL0UT = CELAY(2005.DI LI NT .OUT) 

I F ( TIMfc.Lb'.DELIMT) GO t 0 2 

I F ( ABS(DELOUT-OUT) .LE.tPSl ) GO TO 1 

MAXOUT=0. 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IF(OUT.GT.MAXOUT) maxoijt«out 
DE ROUT- OMEGA *MAXOUT*CCS ( ARG) 

EPS2=0.001*MAX0UT 
EPS3=0MFGA*MAX0UT*CCS( ERG) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 A X I N=C p 0SS(TIME»IN»0.) 

IF( AXIN.NE. 0. ) X I N =A X I N 
AXOUT=CFOSS( TIME .OUT ,0. ) 

IF< AXOUT.NE.O. ) XOUT=AXOUT 
I F ( XIN.GT.XOUT) GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE=?*PI*XTDIFF /PER IOD 
DPHAS2=360.*XT0IFF /PERIOD 

3 APG=OMFGA*TI ME+RPHASE 

I F ( TIME. LE. DILINT) GO TO A 

IF( AB SIDE LOUT-OUT) . LE . E P S2 . AND. ABS ( 0 ILOUT- 3UT . . . 
).LE.EPS2.AND.ABS(DE POUT . . . 

) .LE . EFS3.AND.0UT.GT. 0. . ND . MAX OUT . GT . 0 . ) CALL ENDRUN 
A CONTINUE 
TERMINAL REGION 

PRINT OMEGA , MAXOUT.RPHASc .DPHASE 
CALL PRINT 

END 

PARAM 0MEGA=7. 

CONTRL DELT=0.0A487989 
END 



FIGURE A . 1 

PROGRAM FOR SYSTEM OF EXAMPLE 1 
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PARAM QMEGA=8. 

CONTRL DELT=0. 03926991 
END 

PARAM 0MEGA=9. 

CONTRL DELT=0. 03490659 
END 

PARAM OMEGA= 10 . 

CONTRL DELT=0. 031415927 
END 

PARAM OMEGA= 11. 

CONTRL DELT=0. 02855993 
ENC 

PARAM OMEGA= 12. 

CONTRL DELT=C. 02617994 
END 

PARAM OMEGA= 13 . 

CONTRL DELT=0. 02416609 
END 

PARAM 0MEGA=14. 

CONTRL DELT=0. 02243995 

END 

STOP 



FIGURE A . 1 (CONTINUED) 
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Frequency Response for System of Examples 



EXAMPLE 2 



The system used in this example is identical to the system used 
in the sample program of Figure 3.1. This example shows a second 
method, a variation of that shown in Example 1, whereby the response 
to more than one frequency may be obtained in one job. For this method, 
DELT need not be recalculated for each OMEGA, which is halved each run, 
since halving OMEGA results in doubling the number of iteration inter- 
vals per cycle for a given DELT. To ensure that the criteria for DELT 
(and N) are always met, the highest value for OMEGA for which a response 
is desired in the job is specified in the PARAM card. 

The program is shown in Figure A. 3. The frequency response, 
identical to that obtained for Example 1, is shown in Figure A. 2. 

Since the number of iterations per cycle doubles each run, simu- 
lation time for this method is greater than for the method of Example 
1, where the number of iterations per cycle remains constant for all 
frequencies (so long as the same value for N is used). For a maximum 
allowable simulation time of ten minutes (a local ruling specified for 
the computer on which these examples were run), the responses for a 
maximum of only six frequencies were obtainable on any single job, an 
average of 100 seconds per run. The method of Example 1 is, therefore, 
in general less costly. In general, simulation time may vary consider- 
ably from one system to another and, for any given system, from one 
frequency to another. 
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TITLE LINEAR SYSTEM RESPONSE TO A SINE HAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: 100/S(S+1) 

* FEEDBACK-PATH TRANSFER FUNCTION: -1 
CONST PI=3. 1415927 

PAR AM 0MEGA=40.» LOMEGA= 1 . » A= 1 . » FGA IN= 1. 

CONTRL F I NTI M = 1 OCOO. , DELT=0. 015707964 
INTEG RKSFX 
INITIAL REGION 

0MEGA=0MEGA/2. 

PERI0D=2. *PI /OMEGA 
DELINT=3.*PERI0D 
DILINT=10.*PERI0D 
ERG=PI*(.5-DELT/PERI0D) 

EPS1 =1 • E-C3* A 
EPS2=0. 

EPS3=0 • 

XIN=0. 

XCUT=0, 

XT D I FF = 0 • 

RPHASE=C. 

DPHASE=0. 

MAXOUT=C • 

DEROUT = 10* 

DERIVATIVE REGION 

I N= A*S I NE ( 0. , OMEGA, 0. ) 

ERROR= I N-FGAI N* OUT 

OUT=TRNFR ( 0 . , 2 . , IC , NUM , DEN , ERROR > 

STORAG IC ( 2) , NUM ( 1 I , DEN ( 3 ) 

TABLE IC ( 1-2 I =C* ,0,,NUM( 1 ) = 100. ,DE N ( 1 -3 ) = 1 . , 1 . , 0 . 

DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DEL0UT = DELAY(605, DELINT, OUT I 
D I LOUT = DEL AY (2005, DILINT, OUT ) 

IFITIME.LE.DELINTI GO TO 2 
IF(ABS(DEL OUT- OUT I • LE • EPS 1 1 GO TO 1 
MAXOUT =0 • 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IF (OUT, GT. MAXOUT I MAXOUT=OUT 
PEROUT=OMEGA*MAXOUT*COS( ARG> 

EPS 2=0 ,001 *MAXOUT 
EPS3=0MEGA*MAX OUT *COS ( ERG I 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 AXIN=CROSS (T IME, IN,0. ) 

IF(AXIN,NE,0* ) XI N= AX I N 
AXOUT=CROS S(TIME,OUT,0* ) 

IF (AXOUT.NE.O. ) XOUT=AXOUT 
IF( XIN.GT. XOUT» GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE=2*P I*XTDIFF/PERIOD 
DPHASE=36C,*XTD IFF /PERIOD 

3 ARG=OMEGA*T IME+RPHASE 
IF(TIME. LE. DILINT ) GO TO 4 

IF(ABS( DEL OUT-OUT) . LE. EPS 2. AND. ABS ( DI LOUT-OUT. . . 
).LE.EPS2.AND. ABS ( DEROUT .. . 

).LE.EPS3. AND. OUT , GT. 0. .AND. MAXOUT . GT . 0 . ) CALL ENDRUN 

4 CONTINUE 
TERMINAL REGION 

PRINT OMEGA, MAXOUT, RPHASE , DPHASE 
CALL PRINT 

IF ( OMEGA. GF.LOMEGA) CALL RERUN 

END 

STOP 



FIGURE A. 3 

PROGRAM FOR SYSTEM OF EXAMPLE 2 



30 



EXAMPLE 3 



This example, too, uses the system shown in the sample program of 
Figure 3.1 and in Examples 1 and 2. In this case, however, the system 
is modeled in its open-loop form, rather than the closed-loop form 
employed in the sample and Examples 1 and 2. The program is shown in 
Figure A. 4, and its frequency response, identical to those obtained in 
Examples 1 and 2, is shown in Figure A. 2. 
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TITLE LINEAR SYSTEM RESPONSE TO A SINE WAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: 100/ (S2+S+100 ) 

* FEEDBACK-PATH TRANSFER FUNCTION: ZERO 
CONST PI=3. 1415927 

PARAM 0MEGA=160. ♦LOMEGA*!. »A = 1. 

CONTRL FINTIM=10000. ,DELT=0. 003926989 
INTEG RKSFX 
INITIAL REGION 

0MEGA=0MEGA/2. 

PERI0D=2.*PI /OMEGA 
DELINT=3.*PERI0D 
DILINT=10.*PERI0D 
ERG=P I * ( * 5-DEL T /PERIOD) 

EPS1=1.E-03*A 
EP S2=0 • 

EPS3=0. 

XIN=0. 

X0UT=0. 

XTD IFF=0. 

RPHASE=0. 

DPHASE = 0 . 

MAX0UT=0. 

DER0UT=10. 

DERIVATIVE REGION 

IN=A*$INE(0. , OMEGA, 0. ) 

OUT =TRNFR ( 0. ,2. 1 1 C tNUM ,DEN * I N) 

STORAG IC ( 2 ) « NUM( 1 ) f DEN ( 3 ) 

TABLE IC ( l-2)=0.»0.tNUM(l)=100.»DEN(l-3)=l.»l.»100. 

DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DELOUT=DELAY( 605 .DELINT, OUT) 
DIL0UT=DELAY(2005*DILINT*0UT) 

I F ( TIME. LE. DELINT) GO TO 2 
IF(ABS(DELOUT-OUT).LE.EPS1 ) GO TO I 
MAX0UT=0. 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IF(OUT.GT.MAXOUT) MAXOUT=OUT 
DEROUT=OMEGA*MAXOUT*COS( ARG) 

EPS2=0.00I*MAX0UT 
EPS3=0MEGA*MAX0UT*C0S( ERG) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 AXIN=CR0$S( TIME ,IN,0. ) 

IFUXIN.NE.O.) X I N =A X I N 
AXOUT=CROSS( TIME,OUT f O. ) 

IF( AXOUT.NE.O. ) XOUT=AXOUT 
I F ( XIN.GT.XOUT) GO TO 3 
XTD I FF = XI N-XOUT 
RPHASE=2*PI*XTD IFF /PERIOD 
DPHASE=360.* XTD IFF /PERIOD 

3 ARG=OMEGA*TI ME + R PHASE 

IF( TIME. LE. DILINT) GO TO ^ 

I F ( ABS(DE LOUT- OUT). LE.EPS2. AND. A8S( DILOUT-OUT . . . 

). LE. EPS2. AND. AB SIDE ROUT.. . 

) .LE.EPS3.AND.0UT.GT.0.. AND. MAX OUT . GT. 0 . ) CALL ENDRUN 

4 CONTINUF 
TERMINAL REGION 

PRINT OMEGA » MAX OUT »RPHA$E» DPHASE 
CALL PRINT 

IFIOMEGA.GE. LOMEGA) CALL RERUN 

END 

STOP 



FIGURE A. 4 

PROGRAM FOR SYSTEM OF EXAMPLE 3 
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EXAMPLE 4 



The system used in this example is open-loop, with a transfer 
function of 10/s (s+1) (s+5) . The program is shown in Figure A. 5, and 
the system response is shown in Figure A. 6. 
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TITLE LINEAR SYSTEM RESPONSE TO A SINE WAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: 10/ ( S3+6S2+5S- 10) 

* FEEDBACK-PATH TRANSFER FUNCTION: -1 
CONST P 1=3.1415927 

PAR AM OMFGA=200. » L0MEGA=0. 0 1 , A=l. 

CONTRL FINTIM=10000. , DELT=0. 001 5707964 
INTEG RKSFX 
INITIAL REGION 

0MEGA=0MEGA/2. 

P ER 1 0D=2 »*P I /OMEGA 
DELINT=3.*PERI0D 
DILI NT =10. *PER IOD 
ERG=PI*(.5-DELT/PERI0D) 

EPS1=1.E-03*A 

EPS2=0. 

EPS3=0. 

X I N = 0 • 

X0UT=0. 

XTDIFF=0. 

RPHASE = 0. 

CPHASE=C. 

MAX0UT=0. 

DEROUT =10* 

DERIVATIVE REGION 

I N=A*S I NE ( 0. , OMEGA, 0. ) 

OUT=TRNFR( 0. ,3. ,1 C » NUM, DEN , IN) 

STORAG ICO) ,NUM(1 ),DEN(4) 

TABLE IC( 1-3) =0. ,0. ,0. ,NUM( 1 )=10. , DENI 1-4 )=1 , ,6. ,5. ,0. 
DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DELOUT = DFL AY( 605, DELINT, OUT) 

D ILOUT = DEL AY (2 00 5, DILINT, OUT ) 

IFITIME.LE. DELINT) GO TO 2 
I F ( ABS I DEL OUT- OUT ).LE.EPS1) GO TO 1 
MAXOUT = 0. 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IFIOUT.GT. MAXOUT) MAXOUT=OUT 
DEPOUT =OMEG A* MAXOUT*COS(ARG) 

EPS2=0.0C1*MAX0UT 
EPS3=0MEGA*MAX0UT*C0S( ERG) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 AXIN=C ROSSITIME, I N , 0. ) 

IF! AXIN.NE.O. > XIN=AXIN 
AXOUT=CROSS( TIME, OUT, 0. ) 

I F ( AXOUT.NE.O. ) XOUT=AXOUT 
IF( XIN.GT. XOUT) GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE=2*P I*XTD IFF /PERIOD 
DPHASE=360.*XTD IFF /PER IOD 

3 ARG= OMEGA*T IME+R PHASE 
IFITIME.LE.DILINT) GO TO 4 

I F ( ABSI DEL OUT-OUT J.LE.EPS2. AND. ABSI DI LOUT-OUT. . . 

) .LE.EPS2. AND. ABS ( DEROUT. . . 

J.LE.EPS3. AND. OUT.GT.O. .A ND. MAX OUT • GT • 0 • ) CALL ENDRUN 

4 CONTINUE 
TERMINAL REGION 

PRINT OMEGA, MAXOUT, R PH ASE , DPH AS E 
CALL PRINT 

IFIOMEGA.GE.LOMEGA) CALL RERUN 

END 

STOP 



FIGURE A. 5 

PROGRAM FOR SYSTEM OF EXAMPLE 4 
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Frequency Response for System of Example 



EXAMPLE 5 



The system used in this example is open-loop, with a transfer 
function of (s+1) / (s+10) . To meet DSL/360 criteria for system 
modeling [Ref, 1], the division was carried out, yielding 1 - 9/(s+10). 
The system was then modeled from the diagram shown in Figure A, 7. 




FIGURE A. 7 



Block Diagram of System for Example 5 
The program is shown in Figure A. 8, and the frequency response is 
shown in Figure A. 9. 

The phase of the output in this system leads the input. One 
limitation to the program, discussed in Chapter IV, is that the pro- 
gram computes the supplement of the phase angle rather than the true 
value of the phase whenever the output leads the input, as in this 
example. True phase is used in the plot of Figure A. 9. 
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TITLE LINEAR SYSTEM RESPONSE TO A SINE WAVE INPUT 

* FORWARD-PATH TRANSFER FUNCTION: (S+D/IS + IOI = 

* 1 - 9/IS+10) 

* FEEDBACK-PATH TRANSFER FUNCTION: ZERO 
CONST P 1=3* 1415927 

PARAM OMEGA= 160. ,LOMEGA=l. »A = 1. 

CONTRL F I NT I M=1 0000. ,DELT=0. 003926989 
INTEG RKSFX 
INITIAL REGION 

0MEGA=0MEGA/2. 

PER I0D=2.*PI /OMEGA 
DELINT=3.*PERI0D 
DILINT=10.*PERI0D 
ERG-P I *( . 5-DELT /PERI OD) 

EPS1=1.E-03*A 

EPS2=0. 

EPS3=0. 

XIN = 0 . 

X0UT=0. 

XTD I F F = 0. 

RPHASE =0. 

DPHASE=0. 

MAX0UT=0. 

DER0UT=10. 

DERIVATIVE REGION 

IN=A*SINE ( 0. »OMEGA »0. ) 

OUTA=TRNFR( 0. ,1. , I C , NUM,DE N , I N) 

STORAG IC( 1 ) »NUM( 1 ) ,DEN( 2) 

TABLE IC(1)=0.,NUM(1)=-9.»DEN(1-2)=1.»10. 

OUT=OUTA+IN 
DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DEL0UT=DELAY(605,DELINT,0UT) 

DIL0UT=DELAY(2005,DILINT,0UT) 

IF< TIME. LE. DELINT) GO TO 2 

IF ( ABS( DEL OUT- OUT) . LE . EPS1 I GO TO 1 

MAX0UT=0. 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IFIOUT.GT. MAXOUT) MAXOUT=OUT 
DEROUT=OMEGA*MAXOUT*CGS( ARG) 

EPS2=0.001*MAX0UT 
EPS3=0MEGA*MAX0UT*C0S( ERG) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASF 

2 AXIN=CROSS(TIME ,IN,0. ) 

IF( AXIN.NE.O.) XI N=AXI N 
AX0UT=CR0SS(TIME,0UT,0.) 

IF( AXOUT.NE.O. ) XOUT=AXOUT 
IFI XIN.GT.XOUT) GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE-2*PI*XTDIFF /PERIOD 
DP HA SE=360.* XTD IFF /PERIOD 

3 ARG=OMEGA*TI ME+RPHASE 
IFITIME.LE. DILINT) GO TO A 

IF ( ABS( DE LOUT-OUT) .LE.EPS2.AND. ABS ( D I LOUT-OUT .. . 

) .LE.E PS2. AND. AB S ( DE ROUT.. . 

) .LE.EPS3.AND.0UT.GT.0. . AND. MAXOUT . GT. 0 . ) CALL ENDRUN 
A CONTINUE 
TERMINAL REGION 

OMEGA, MAXOUT , RPHASE ,DPHASE 
CALL PRINT 

IF ( OMEGA .GE . LOMEGA ) CALL RERUN 



PRINT 



END 

STOP 



FIGURE A. 8 

PROGRAM FOR SYSTEM OF EXAMPLE 5 
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FIGURE A. 



EXAMPLE 6 



The system modeled in this example is a magnetic tape transport 
mechanism which employs a vacuum column for tape loop control. The 
equations representing the system have been linearized for use with 
the program. 

The input is an AC signal to a DC motor, and the output is the 
velocity of the tape before it passes over the capstan. 

A diagram of the system is shown in Figure A. 10, and the program 
containing the system modeling equations is shown in Figure A. 11. 

The frequency response of the system is shown in Figure A. 12. 

It will be noted that a resonance peak occurs at a frequency of 
approximately 20,000 radians per second. All attempts to determine 
frequency response at or near this frequency were unsuccessful because 
the steady-state criteria were not met within ten minutes, the maximum 
simulation time available locally on the computer on which these 
examples were run. 
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FIGURE A. 10 

Diagram of System for Example 6 
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TITLE LINEAR SYSTEM RESPONSE TO A SINE WAVE INPUT 
* SYSTEM- -LINEARIZED C YBE RNE T CAPSTAN MOTOR , TAPE , AN 

CONST PI=3. 1415927,... 

MA=19. 245F-06, K4= 1 1 9. 8 , C4=8 . 5 E- 03 , RCAP= . 75 , . . . 
LM=20.0E-06,RMOTOR=0.6,KBEMF=. 0397, M0=2. 1 15E-06,.. . 
K3= 13 5. ,C 3=. 0092, K0= 1 091., C0=. 0255, M3=14.15E-06, ... 

J 2=22. 5E-C6,K1 2 = 6000. ,C 12=0. 008 , J1 = 50. E-06 , . . . 

KT=0. 353, C 1=14. F-06,. .. 
FM=C.0»F0=0.0,F3=0.0»F4=0.0,T0RK=0.0 
INCON ICM=C. , ICC 1=0. , IC02=l.E-30, IC11=0. , I C 12 = 1 . E-30 , . . . 

IC22=C. , IC31=0. , IC32=l.E-30, IC41=0. , I C42= 1 . E- 30 , . . . 
IC21=0. 

PARAM 0MEGA=450C0. ,L0MEGA=1. ,A=1. 

CCNTRL F I NT I M=10000. , DELT=0 .000 013^626 
INTEG RKSFX 
INITIAL REGION 

OMFGA=OMEG A/2. 

PERI 00=2. *PI /OMEGA 
DELINT=3.*PERI00 
DI L I NT = 1 0. *P ER I OD 
ERG=PI*(.5-DELT/PERI0D) 

EPS1=1 . E-03*A 
EPS2=0. 

EPS3=0. 

X I N=C. 

XOUT =0. 

XTD I FF= 0 . 

RPHASE=0. 

DPHASE=C. 

MAXOUT=0. 

DEROUT = 10. 

DERIVATIVE REGION 

IN=A*SINF(0. .OMEGA, 0 . ) 

VMOT 0R= I N 

LM I DOT = VMO TOR-RMO TOR* IM-K8EMF*TH1D0T 
I MDCT= ( l./LM)*LMIDOT 
I M=I NTGRH IC M , I MOOT ) 

MCTORQ=KT* I M-C 1*TH 1D0T- SI GN ( FM , TH1D0T ) -K1 2*. . . 

(TH1-TH2 )-C12*(THl DOT-T H2D0T ) 

TH ID 2= I 1./ J1 )*MOTORQ 
TH1D0T= INTGRLI IC12.TH102) 

TH1= I NTGPL ( I Cl 1 ,TH1 DOT ) 

SUM4=T0RK- T2-S I GN( F4, X4D0T ) 

X4D2=( l./M4)*SUM4 
X4D0T= INTGRLI IC42,X4D2) 

X4= I NT GRL( IC41.X4D0T) 

SUMX42 = X4- X2 

D0TX42 = X4D 0T-X2D0T 

T 2=C 4*00TX 42+K4*SUM X42 

X2D0T=X1D0T 

X2=X1 

T12=T1-T 2 

CAPT0R = K12* (TH1-TH2 I+C12* I T HI DOT- TH 2 DOT )-RCAP*T12 
TH2D2= ( 1. / J2 ) *C APTOR 
TH2D0T= INTGRLI IC22.TH2D2) 

TH2= I NTGRL ( IC21 ,TH2D0T ) 

T 1=K0*SUMX01+C0*D0TX01 
SUMX01 = XI- XO 
D0TX01 = X1DOT-XODOT 
X1D0T=RCAP*TH2D0T 

X I = RCA P*TH2 

SUMO=K 3* X3 + C3* X3D0T + K0* XI +C0*X1 OOT- ( C0 + C3 ) *XODOT-. • • 

I K0+K3 )*XC-SIGN(FQ» XODOT) 



FIGURE A. 11 

PROGRAM FOR SYSTEM OF EXAMPLE 6 
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X0D2= ( 1 • /HO ) *SUMO 
XODOT= INTGRLUC02, XOD2) 

XO = I NTGR L( ICOl .XODOT) 

SUM3=K3*(XO-X3)+C3*< XODOT-X3DOT ) -TORK-S IGN ( F3 , X3DOT ) 

X3D2=(1./M3)*SUM3 

X3DOT=INTGRL< IC32,X3D2) 

X3= INTGRL( IC31 » X3DOT ) 

OUT=XODOT 
DYNAMIC REGION 

* DETERMINATION OF STEADY-STATE CONDITION 

DE LOUT =DELAY(605,DELINT .OUT ) 

DIL0UT=DELAY(2005»DILINT, OUT ) 

IF(TIME.LE.DEL INT) GO TO 2 

I F ( ABS ( DEL OUT-OUT ) • LE • EPS 1 ) GO TO I 

MAX0UT=0. 

GO TO 2 

* DETERMINATION OF OUTPUT STEADY-STATE PEAK MAGNITUDE 

1 IF(OUT.GT.MAXOUT) MAXOUT=OUT 
DEROUT =OMEGA*MAXOUT*COS ( ARG) 

EPS2 =0*001* MAX OUT 
EPS3=0MEGA*MAX0UT*C0S( ERG ) 

* DETERMINATION OF OUTPUT STEADY-STATE PHASE 

2 AXI N=CROSS (TIME, IN, 0. I 
IF(AXIN.NE.O.) XIN=AXIN 
AXOUT=C ROSS ( T I ME »OUT »0* ) 

I F ( A XOUT.N E • 0* ) XOUT=AXOUT 
IFIXIN.GT.XOUT ) GO TO 3 
XTDIFF=XIN-XOUT 
RPHASE=2*P I *XTD IFF /PERIOD 
DPHASE=360**XTD IFF /PERIOD 

3 ARG= OMEGA* T I ME + RPHASE 
IFITIME.LE. DILINT) GO TO 4 

I FIABSI DEL OUT-OUT ) *LE*EPS2. AND. ABS ( DI LOUT-OUT*. . 

) • LE.EPS2* AND. ABS( DEROUT. .. 

).LE.EPS3.AND.0UT.GT.0..AND.MAX0UT.GT.0. ) CALL ENDRUN 

4 CONTINUE 
TERMINAL REGION 

PRINT OMEGA, XO, OUT, MAXOUT, XI, XI DOT, X4, X4D0T ,RPHASE, DPHASE 
CALL PRINT 

IF(OMEGA.GE.LOMEGA) CALL RERUN 

END 

STOP 



FIGURE A. 11 (CONTINUED) 
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Frequency Response for System of Example 



APPENDIX B 



A program for generating a table of DELT versus OMEGA for an N of 
twenty is shown in Figure B.l. To generate a table for any other value 
of N, the user need only specify the desired value of N in the program. 
The program may also be modified to obtain values for DELT for intervals 
of OMEGA and a maximum value of OMEGA other than the one radian per 
second and the 10,000 radians per second, respectively, shown in the 
program. 



C GENERATION OF TABLE OF DELT VERSUS OMEGA FOR N=20 
C 

WRITE (6,1) 

1 FORMAT( IX, 'TABLE OF DELT VERSUS OMEGA FOR N=20',//,1X, 
*' OMEGA' ,9X, 'DELT' ,/) 

PI=3. 1415927 
OMEGA-1 . 

AN=20 . 

DO 3 1=1,10000 
DE LT=2 . *PI/ (AN*OMEGA) 

WRITE (6, 2) OMEGA, DELT 

2 FORMAT (F7 . 0,E 18. 7) 

0MEGA=0MEGA+1 . 

3 CONTINUE 
STOP 
END 



FIGURE B.l 



Program for Generating Table of DELT versus OMEGA 
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